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FIELD OF THE INVENTION 

The present invention generally relates to aviation electronics, and more particularly 
relates to multi-sensor navigation systems, and even more particularly relates to methods 
and systems for fault detection and fault exclusion in multi-sensor navigation systems. 
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BACKGROUND OF THE INVENTION 

In the past, multi-sensor navigation systems have been used with much success. 
These multi-sensors typically use global positioning system (GPS) receivers in combination 
with other inertial sensors to generate position information. However, due to the 
inadequate signal-in-space integrity currently provided by GPS for civil aviation, it is 
generally necessary to also employ real-time fault detection and exclusion (FDE) schemes, 
especially when the multi-sensor system is used for primary means of navigation. Since 
the current assumptions about the GPS constellation and its measurement accuracy often 
do not provide high enough availability of FDE to meet the primary means objectives for all 
phases of flight, integration of GPS with inertial sensors for integrity monitoring has 
received considerable attention. See Brenner, M., "Integrated GPS/lnertial Fault Detection 
Availability", ION GPS-95. Palm Springs, CA, Sept. 2-15, 1995, pp. 1949-1958, hereafter 
referred to as "Ref [1]." Also see Vanderwerf, K., "FDE Using Multiple Integrated 
GPS/lnertial Kalman Filters in the Presence of Temporally and Spatially Correlated 
Ionospheric Errors", ION GPS-2001, Salt Lake City, UT, Sept. 11-14, 2001, pp. 2676-2685, 
hereafter referred to as "Ref [2]". Also see Diesel, J., Luu, S., "GPS/IRS AIME: Calculation 
of Thresholds and Protection Radius Using Chi-Square Methods", ION GPS-95, Palm 
Springs, CA, Sept. 12-15, 1995, pp. 1959-1964. hereafter referred to as "Ref [3]." Prior art 
GPS integrity algorithms can be classified into two broad categories based on the method 
employed for computing the navigational solution ~ snapshot approaches and filtered 
approaches. Snapshot approaches are generally based on least-squares (LS) methods, 
while filtered approaches generally utilize multiple Kalman filters with different fault 
hypothesis models. 
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In addition to the snapshot vs. filtered class'fication, the prior art GPS integrity 
algorithms can also be classified into two different categories based on the characteristics 
of their test statistics for the FDE function -- range domain methods vs. position domain 
methods. For example, the solution separation algorithms described in Refs. [1], [2], and 
in Bnjmbacl<, B. D., Srinath, M. D., "A Chi-Square Test for Fault Detection in Kalman 
Filters", IEEE Transactions on Automatic Control, Vol. AC-32, No. 6, June, 1987, pp. 552- 
554, hereafter referred to as "Ref. [4]" are position domain methods, and the measurement 
residual algorithm in Ref. [3] is a range domain method. 

The prior art multi-sensor systems and FDE algorithms have several common 
characteristics and requirements. The FDE function is typically required to ensure the 
integrity of the navigation solution and prevent the use of hazardous and misleading 
information. The prior art FDE function typically consists of two distinct parts ~ fault 
detection and fault exclusion. The purpose of prior art fault detection has been to detect 
the presence of an unacceptable error, and the purpose of the typical prior art fault 
exclusion has been to identify and exclude the culprit causing the unacceptable error with 
high confidence. The FDE function is also often required to provide a statistical bound, 
horizontal protection level (HPL), which the FDE function often guarantees that the 
horizontal position error can only exceed it prior to fault detection within the specified 
probability of missed detection (Pmd)- 

After a fault is detected in prior art systems, the HPL is often operationally of little 
value because HPL no longer bounds the true en-or at the specified statistical level. The 
horizontal uncertainty level (HUL) can provide a bound on the horizontal position error after 
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fault detection and before the correct fault exclusion can be performed. Thus, it typically 
enables the pilots or airborne equipment to determine whether the navigation solution is 
still acceptable at a given phase of flight. The operational differences between HPL and 
HUL are discussed in Young, R. S. Y., McGraw, G. A., Driscoil, B. T.. "Investigation and 
Comparison of Horizontal Protection Level and Horizontal Uncertainty Level in FDE 
Algorithms", ION GPS-96, Kansas City, MO, SepL 17-20, 1996, pp. 1607-1614, hereafter 
referred to as "Ref. [5]. " 

While these prior art multi-sensor systems and FDE algorithms have enjoyed 
success in the past, they have some drawbacks in certain circumstances. Often these 
systems fail to provide a sufficiently precise probability level for the detection threshold and 
HPL. Additionally, in some prior art systems, it is necessary to make real GPS 
measurements to achieve a requisite level of predictability. 

Consequently, there exists a need for improvement in the FDE algorithms for multi- 
sensor navigation systems. 
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SUMMARY OF THE INVENTION 

It is an object of the present invention to provide a multi-sensor GPS navigation 
system which provides data which is sufficient to serve as a primary navigation source for 
commercial air transport aircraft. 

It is a feature of the present invention to utilize a fault detection algorithm which is 
based on normalized solution separation which uses a normalized approach to determine 
statistical properties of a test statistic. 

It is an advantage of the present invention to provide enhanced ability to obtain 
sufficient probability levels for the detection threshold and HPL. 

It is another advantage of the present invention to predict HPL without the need for 
making real GPS measurements. 

It is another object of the present invention to provide a tightly integrated 
GPS/inertial sensor system which can be used for the primary means of navigation from 
oceanic/remote down to non-precision approach phases of flight. 

It is another feature of the present invention to employ a hybrid algorithm which 
combines two independent and parallel methods for fault detection and fault exclusion. 

It is an advantage of the present invention to provide for speedier determination of 
fault exclusions in some circumstances. 

It is another feature of the present invention to provide an indication of HUL for a 
combined GPS/inertial navigation system. 
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It is another advantage to provide additional valuable information to the pilot. 

It is yet another feature of the present invention to provide for parallel independent 
processing schemes for fault exclusion, which schemes include a post update residual 
monitoring scheme, together with a well-known least squares scheme. 

It is yet another advantage of the present invention to eliminate the need for making 
at least six GPS measurements to exclude a fault in circumstances having medium to fast 
failure rates. 

The present invention is a system, method and apparatus designed to achieve the 
above-mentioned objects, include the earlier-listed features and provide the already 
articulated advantages. 

Accordingly, the present invention is a combined GPS/inertial sensor system which 
uses independent fault detection and fault exclusion schemes, as well as improved 
methods for fault detection by determining B(k) (the covariance of solution separation 
vector) in a more efficient way and a post-update residual monitoring scheme for fault 
exclusion. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

Figure 1 is a block diagram of the GPS/INS integration architecture of the present 
invention. 

Figure 2 is the processing diagram of the FDE algorithm for the integration 
architecture of Figure 1 . 

Figure 3 Is a conceptual illustration of Eq. (25) for ||m||»0. 

Figure 4 is a conceptual illustration of Pbias for the 2-DOF case. 

Figure 5 is a flow diagram of the fault exclusion function where the blocks with bold 
solid lines are designed to meet the first requirement for fault exclusion, and the block with 
bold dashed lines is designed to satisfy the second requirement. 
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DETAILED DESCRIPTION 

Now referring to the Figures, where like numerals refer to like matter and text 
throughout, and more particularly referring to Figure 1, generally designated 100, there is 
shown a block diagram of the GPS/INS integration architecture of the present invention. 
This is an open loop integration such that there is no feedback to the INS 106 or the GPS 
receiver 104 from the integration processor 110. Air data 102 is shown to be provided to 
both the INS 106 and the integration processor 1 10 which outputs a corrected position 112, 
as well as other data including HPL, HUL, shown in Figure 2. The proposed Kalman filter 
is executed in the integration processor 110. Also in Figure 1, 108 represents an altitude 
line; 1 14 represents a pseudoranges and SV ECEF positions line; and 116 represents a 
position, velocity and attitude line. 

Now referring to Figure 2, there is shown the architecture, generally designated 200, 
of the integration processor 110. There are multiple navigation solutions for the filtered 
FDE algorithm. The filtered FDE algorithm uses one Full-Filter 202, multiple Sub-Filters 

203 and 204 and one least-squares navigation method 206. The Full-Filter 202 and the 
least-squares method 206 process all available measurements, while the Sub-Filters 203- 

204 each excludes one of available measurements. The fault detection function 208 uses 
information from the Full-Filter 202 and all Sub-Filters 203-204. The fault exclusion 
function 210 uses information from available measurements, the least-squares method 
206, and all Sub-Filters 203-204. Note: is the estimated horizontal position error vector. 
Also in Figure 2, 201 represents an INS position line; 205 represents a re-initialization line; 
and 207 represents a corrected position line. 
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The present invention can be better understood by first understanding and then 
contrasting it with the prior art GPS integrity algorithms, where test statistics for both the 
fault detection function and the fault exclusion function are in the same domain; i.e., both in 
the position domain or both in the range domain. In such prior art systems, the processing 
of the fault exclusion function is usually started after a fault is detected by the fault 
detection function. 

In contrast, for the FDE algorithm of the present invention, the test statistics for the 
fault detection function 208 and the fault exclusion function 210 are in different domains; 
i.e., the former is in the position domain, and the latter is in the range domain. Therefore, 
for some failure types, the test statistic for the fault exclusion function 210 may detect the 
presence of the failure earlier than that for the fault detection function 208. To take 
advantage of this characteristic and achieve better FDE performance, it is proposed to 
treat the fault detection function 208 and the fault exclusion function 210 as two 
independent processes that run in parallel. 

Kalman Filter Structure 

The simulation results presented herein are based on design with 31 Kalman error 
states for each filter. The states are as follows: 

5xo.j : Angular position errors 

: Altitude error 
SX3.5 : Velocity errors 
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gdc^g : Attitude and heading errors 

8x9.,! : GPS receiver clock bias, clock frequency error and clock frequency rate error 

: Baro altitude bias 
Sxi3.,5 : Accelerometer biases 
5*16-18 ■ Gyf'o biases 
2*19-30- GPS pseudorange biases 

The first 9 error states are referenced to the local navigation frame. The 
accelerometer and gyro biases are referenced to the body frame. The GPS pseudorange 
biases are along the Line-Of-Sight (LOS) of each satellite in view. The accelerometer and 
gyro biases are each modeled as a first-order Markov process. The GPS pseudorange 
with Selective Availability absent is also modeled as a first-order Markov process. The 
pseudorange measurement noise was modeled as a zero mean, Gaussian white 
sequence. To balance between the computational burden and navigational accuracy, the 
measurement update period Tb for all Kalman filters was chosen to be greater than one 
second. 

FAULT DETECTION PROCESSING 

The normalized solution separation algorithm for the fault detection function was 
adapted from the chi-square test developed by Brumback and Srinath in Ref. [4]. The 
assumptions, processing, and essential properties of the algorithm are discussed below. 
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System and Filter Models 

Assume the following discrete state space model for the true system error state x(k) 
and the measurement z(k) for the Full-Filter. 

ji:(k) = 0(k,k-lMk-l)+<D(k) 

z(k) = ^,(k)x(k) + Uo(k)+*(k)u(k,ko) ^ ^ 

where x(k) Is a £xl vector and z(k) is a nxl vector. a> (k) and Oq (k) are 
independent, zero means, Gaussian white sequences having covariance of Q(k) and 
i!o(k) , respectively (the subscript "0" denotes parameters associated with the Full-Filter). 

The initial state x(0) is a Gaussian random vector independent of g) (k) and (k) and has 
mean and covariance . The nxl b(k) vector models the failed event, i.e., 
b(k) = b(k>j . u(k,ko) is the unit step function at k^ , i.e., u=l for k > ko , where k^ is the 
time at which the failure occurs. 

When the jth measurement is excluded, the measurement vector in Eq. (1) is 
modified as 

z,{k) = E,zik) (2) 
where E-=I^ -e-ej and =[0,...,0,l,0,...,Of , which is the jth unit standard basis in 

The initial condition and covariance for both Full-Filter and jth Sub-Filter are 
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Xo(0/0) = x.(0/0) = E{x(0)} = x, (3) 

Po(0/0) = P.(0/0) = Po (4) 

Since the jth Sub-Filter uses one less measurement than the Full-Filter, there is one 
pseudorange bias error state used by the Full-Filter, but not by the jth Sub-Filter. However, 
with Eq. (3), this unused pseudorange bias error state in the jth Sub-Filter is initialized to be 
equal to the Full Filter's counterpart. It can be shown that this approach has no effect on 
the performance of the jth Sub-Filter since this SV is not used. However, for real-time 
software, this unused bias state can be dropped for the jth Sub-Filter to save computational 
time. 

Define the estimation errors, and , as: 

xi{k)=xi{kncyx\k) (5) 

x^{k) = x^{\dkyx\k) (6) 

where xl is the horizontal position error estimate (a 2x1 vector) from the Full-Filter, 
x]* is the horizontal position error estimate (a 2x1 vector) from the jth Sub-Filter and x^ is 
the true horizontal position error. 

Derivation of Nornnalized Solution Separation Test Statistics 

The solution separation vector between the Full-Filter and the jth Sub-Filter is 
defined as 
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fi, (k) = ~xl (k) - x] (k) = xl (k/k) - x] (k/k) (7) 

Prior to the onset of the failure, I.e., k < k^ , E|^j(k)| = 0 since each filter is linear 
and each estimate is unbiased. The covariance of the solution separation vector is 

If j (k) = Po" (k/k) - (k/k) - PjS (k/k) + Pi (k/k) (8) 
where 

/>-(k/k) = E{iJ(io^r| 

i^.^(k/k) = E{if(if)'} (9) 
^S(k/k) = E{xf(.i„^)^}=(i>;(k/k))^ 

However, it is shown in a dedicated section below tliat in fact (k/k) = (k/k). 
Therefore, 

(k) = Pf (k/k) - P^ (k/k) (1 0) 

This simplification yields considerable computational savings as compared to the 
dual covariance propagator presented in Refs. [1] and [2]. 

The normalized solution separation test statistic is given by 

X-(k) = /^^(k)5;(k))ff.(k) (11) 

where denotes the Moore-Penrose generalized inverse, as discussed in Horn, 

R. A., Johnson, C. R., Matrix Analysis, Cambridge University Press, 1985, hereafter 

-13- 
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referred to as Ref [6] of B- and X-(k) is chi-square distributed. A fault is declared if X^l^ik) 

> TD2 for any j, where TD2 is the detection threshold. To be consistent with the detection 
criterion used in Brown, R. G., Chin, G. Y., "GPS RAIM Calculation of Thresholds and 
Protection Radius Using Chi-Square Methods— A Geometric Approach", Red Book Series, 
Vol. V, of the Institute of Navigation on Global Positioning System, The Institute of 
Navigation, Alexandria, VA, 1998, pp. 155-178, hereafter referred to as Ref [7]. ^^(X) 
instead of X^(k) is chosen as the test statistic in this paper. Since the relationship between 
?iy^(k) and X^(k) is monotonic, the fault detection algorithm can be developed based on 
the chi-square statistics and the appropriate conversion is used to obtain the corresponding 
result for xy^(k) . The determination of the TD2 will be discussed below. 

Normalized Solution Separation Detection Threshold: Full-Rank Case 

In this case, the Moore-Penrose inverse in Eq. (1 1) becomes a regular matrix 
inverse. However, for reasons of numerical stability in on-line computations and 
consistency with the general case of B^{k) being rank deficient, the computation in Eq. 
(1 1 ) is implemented as follows: 

Denote the eigenvalue decomposition of J?, (k) as 

5,(k) = F3(k)5.(k)F/(k) (12) 
with 
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■<^.,Bj(k) 0 

0 <'2.Bi(k) 



(13) 



F3(k) = [v,(k) v,(k)] (14) 

where v, (k) and (k) are the eigenvectors of (k) and cr,3j > a^^ . The 
normalized solution separation vector is given by 

j.,(k) = 5:'/^(k)F/(k)^.(k) (15) 
with 



5:-(k) = 



0 



(16) 



and the detection statistic becomes 



^(k) = J^;(k)j^,(k) 



(17) 



In the full-rank case, (k) is chi-square distributed with two degrees of freedom 

and the selection of TD2 to meet the required false alarm probability is straightfonA/ard. 
Because there are multiple filters and any of the filters can cause a false alert, TD2 needs 
to be determined from Pfa/N, where Pfa is the specified false alert rate and N is the 
number of available measurements. 

Prior solution separation FDE approaches in the literature (e.g., Refs [1] and [2]) 
have used the unnormalized test statistic jfi. (k)|| along with a one-dimensional Gaussian 
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statistics and standard deviation of ai,Bj to compute the detection threshold and HPL. 
Since is not a Gaussian variable, the use of the 1-DOF Gaussian statistics in 

determining the detection threshold and HPL may result in an underestimate of the actual 
error. 

Normalized Solution Separation Detection Threshold: Rank-Deficient Cas 

When i?j(k) is either rank-1 or ill-conditioned, the inverse of 5j(k) does not exist or 

its computation will lead to numerical difficulties. This numerical stability issue could be 
resolved if the detection statistic is treated as 1-DOF Gaussian distributed (i.e., ignore enior 
contribution along the (k) direction) by computing (k) through the Moore-Penrose 
generalized inverse: 

5:-(k) = 

However, this approach compromises false alert and missed alert probabilities for 
the ill-conditioned case. This is because, unlike the least-squares case where both the 
solution bias and noise lie only in the (k) subspace, in the filtered case there are noise 

and parts of the solution separation bias that are along the (k) direction, which would 

not be accounted for if the detection statistic is treated as 1-DOF Gaussian distributed. 
Therefore, it is proposed to modify Eq. (16) and detection threshold as follows to both 
preserve numerical stability and provide better false alert and missed alert protection than 
just using 1-DOF Gaussian approximation. 



0 
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In the ill-conditioned or rank-1 case, it is proposed to use 

5:-(k). 

A fault is declared if 

;.f(k)>(Ka,3.(k)/a,3.(k)y''TD,, foranyj (20) 

where TDj is determined from 1-DOF Gaussian statistics to meet the probability of 
false alarm Pfa/N. When J?j(k) is a rank-1 matrix, there is no noise in the V2(k) subspace. 
Therefore, assigning the second diagonal term in Eq. (19) as o~^^^(k) does not introduce 
additional noise on X.(k). In addition, Eq. (20) provides the exact false alert probability. 
When B.(k) is rank-1 or ill-conditioned, the second diagonal term in Eq. (19) enables the 
component of the solution separation bias vector along the (k) direction to be properly 
scaled in Xj(k) so it can be accounted for in the computation of the estimated horizontal 
position error later through the multiplication of <y}3j(k) in Eq. (34). Therefore, with Eq. 
(19), the exact magnitude of the solution separation bias vector is captured in X-(k) and 
accounted for in the estimated horizontal position error, which the 1-DOF approximation 
does not since it ignores any error contribution along the V2(k) direction. When B-(k) is ill- 
conditioned, with TDj scaled as in Eq. (20), a portion of the noise in the V2(k) direction is 
accounted for, which the 1-DOF approximation does not. 



0 
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HORIZONTAL PROTECTION LEVEL 

If the horizontal position error exceeds the HPL, it can be due to a measurement 
fault incorporated into the navigational solution or the rare normal fault-free event. See 
"Minimum Operation Performance Standard for Global Positioning SystemAA/ide Area 
Augmentation System Airborne Equipment", Appendix R, RTCA/DO-229C, hereafter 
referred to as Ref. [8]. Therefore, HPL is determined as follows: 

HPL = max {HPLho, HPLh, } (21 ) 

where HPLjjo is based on the rare normal fault-free hypothesis and HPLm is based 
on the fault-in-progress hypothesis. The algorithms for HPL^ and HPLj,o are derived 
below. 

HPL Based on the Fault-ln-Progress Hypothesis (HPLhi) 

The estimated horizontal position error vector from the Full-Filter at the fault 
detection can be expressed as follows: 

^o' (kio) = (I^td) + bias, (kTD)+ ©0 (k^j,) (22) 

where k^^^ is the time at fault detection, ©oCk^i^) is the noise vector with zero mean 
from the Full-Filter and bias^i^-^^) is the horizontal position bias vector of the Full-Filter. 

Since the solution from the Full-Filter is output for use by external airborne 
equipment, the true horizontal position error vector is defined as follows: 



"18- 



PATENT 

Attorney Docket No. 02CR109/KE 

Err{k,^) = {k^^yx' {k^^) = biaSo (k,o) + <0o (k^^) (23) 

Err{kjj^) cannot be computed directly because both biasQ(k^j^) and cDoCk^o) are 
unknown. Therefore, we need to estimate Err(kj^) based on some parameter that can be 
computed directly, namely the solution separation vector, in Eq. (7). 

Before proceeding to the HPLhi algorithm, we introduce the following notation to 
facilitate the mathematical derivation. 

^(®>Pmd) 'S the radius such that Prob{||(D|| > i?((D,Pj^i,)}=Pj^j, (24) 
i?(<o,iii,PMD) is the distance such that Prob{||iii+(D||>||iii||+i?((D,iii,Pj^)}=Pj^ 

(25) 

Note that 7?(o, 0.5) is the standard circular error probability for a two-dimensional 
random variable co. When ||iw||=0, R{ai,P^)=R{G},m,?^). 

Now referring to Figure 3, generally designated 300, 301 represents a portion of a 
noise cloud with Pmd probability; 303 represents R(a),m,PMD); and 305 represents a (1-Pmd) 
error bound of E. Also in Figure 3, an illustration of Eq. (25) for \\m\\ » 0 is shown. With 

Eqs. (23) and (25), the ideal horizontal position error bound, a, can be written as 

a = \\bias,ik^)\\^R{(s>,(k^^)M^^^ (26) 
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Since biaso(kja) is unknown, the solution separation vector fi. is used to estimate It 

with (1-Pmd) probability. At the fault detection, the solution separation vector for the jth 
Sub-Filter not using the failed measurement can be written as 

fi^=xl -if = Woso +©0 -fflj =bias, -®; (27) 

where a. Is the noise vector from the jth Sub-Filter and tn'. is the component of <o. 
independent of (d,, (Note: the index is suppressed here and below for notational 
convenience). 

As shown in Eq. (27), using fi^ to estimate biaSf, introduces an additional noise term 
cDj. that does not exist in Eq. (26). To properly account for the effect of this additional noise 
term, Eq. (26) is re-arranged as 

a = \\bias,\\+R(ia]Mas,,'P^) + R{&^,bias,,f^^)-R^ (28) 

Based on Eq. (28), the algorithm for the ideal horizontal position error bound a is 
divided into two parts. The first part is to estimate the first two terms at the right hand side 
of Eq. (28) due to bias^ and uncorrelated noise coj with fi. . The second part is to estimate 

the last two terms at the right hand side of Eq. (28) due to the noise cOq and Oj . 

Horizontal Position Error Due to Bias and Uncorrelated Noise 

The mean of the normalized solution separation vector and the test statistic can be 
written as: 
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(29) 



(30) 



where the noise components, yj (1) and yj(2) , are independent with N(0, 1) and the 
means for yj.(l) and yj(2) are m^^ and lUj^, respectively. 



distributed. Therefore, we can compute the noncentrality parameter that makes the noise 
inside the detection threshold to be Pmd- With normalization, this noncentrality parameter 
is the "Pbias" in parity space FDE algorithms Ref. [7]. 

Now referring to Figure 4, generally designated 400, 401 represents a noise cloud 
with Pmd inside detection threshold; 402 represents Pbias; and 403 represents a 
normalized solution separation (y) domain line. There is shown an illustration of Pbias for 
the 2-DOF case. So, for any bias vector with a magnitude greater than Pbias, it will be 
detected with (1-Pmd) probability. 404 represents the detection threshold. 

The relationship between Pbias and its corresponding bias vector, , in the 
horizontal position domain can be defined as: 



With an orthonormal rotation, X- can be shown to be noncentral chi-square 



Pbias 



(31) 



Pbias = ,Jl^ 



(32) 
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After determining the Pbias in the normalized solution separation (y) domain, we 
need to convert it back to determine ||/pbias|| since our main interest is the horizontal 



position error, /"pbias be computed as 



llnbiJhV^f^ (33) 

To compute Eq. (33), we need to know the Y. vector defined In Eq. (31 ). Since the 
exact orientation of the ^Vbias vector is unknown, we cannot determine Y^. Nonetheless, we 
can compute the upper bound of Eq. (33) as follows: 







0 




0 


^2,Bj_ 



l^.<Pbias7^ =HPE_B. (34) 



where a, g, is the largest eigenvalue of the Sj matrix and HPE_Bj is the horizontal 
position error bound of ||i^pbias|| with a specified probability of missed detection (Pmd). 



Since the worst case assumption is used in Eq. (34) and p^^bias^ -Oj, it can be 

shown 

HPE_Bj >||Aia5j| + i?(cD;,WaSo,PMD) (35) 
Horizontal Position Error Due to Noise Only 

The last two terms at the right hand side of Eq. (28) can be related to the radii 
determined with the circular error probability calculation as 



PATENT 

Attorney Docket No. 02CR109/KE 

R{<a„bias„F^) = R{(0„V^)-e, (36) 

R ((d;, W«5„,P^ ) = («>'.,P^ ) - 83 (37) 
It can be shown that e, > 0 and > 0 . 

Now, by substituting Eqs. (35)-(37) into Eq. (28), Eq. (28) can be re-written as 

a<HPE_Bj + /?(<Do,PMo)-i?(<D;,PMD)+S2-e, (38) 

R{a)o,f^o) in ^Q- (38) can be decomposed as 

^KPmd) = ^K.Pmd)-83 (39) 
Since Gi.=tdo + (a'., S3 >0. With Eq. (39), Eq. (38) can be re-arranged as 

a<HPE_Bj+i?((Dj,P^)-i?((D;,P^i,) + s,-s,-S3 (40) 

Now, let us examine the error term, -s, -83 , in Eq. (40). For the limiting case 
Ejoj = 0 , 82 = 0, 83 = 0, and 8, ^ 0. Therefore, for the limiting case, a will be bounded 
by 

a<HPE_B,. + i?((D..P^)-i?(fl);,P^) (41) 

It has been found by simulation that Eq. (41) appears to be always valid. Therefore, 
Eq. (41) provides the required statistical bound for horizontal position errors. 
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The last two terms in Eq. (41) can be computed through their respective noise 
covariance and B- by the following two steps. 

(1) The noise on the X and Y axes in the horizontal navigational coordinates are 
likely to be correlated, i.e., the covariance has non-zero off-diagonal terms. Therefore, 
a coordinate transformation should be used to transform the X-Y plane to a principal 
X1-Y1 plane such that noise components in X1 and Y1 axes are not correlated. 

(2) After the coordinate transformation, the noise components in X1 and Y1 axes 
are independent Gaussian noise with, in general, unequal standard deviations. The 
upper bound of the horizontal position error for a specified probability can be computed 
by using circular error probability (CEP) table look-up, Beyer, W. H., "Handbook of 
Tables for Probability and Statistics", The Chemical Rubber Company, Cleveland, Ohio, 
1966, pp. 146-148 hereafter referred to as Ref. [9] or a polynomial approximation. For 
our applications, a CEP table look-up suffices since a constant Pmd is required. 

Following the above two steps, ^^((Oj,?,^) can be computed as 



where j and a^ - are the eigenvalues of P- , a, j > a2 j , and Kcepj is a function of (1- 




(42) 



Pmd) probability and [^x/J^2,i)^^ CEP look-up table. 



^(©j^PMo) is computed as 



PATENT 

Attorney Docket No. 02CR109/KE 
HPE.NBj =i?(a>;,P^)= Kcep,Bi7^ (43) 

where c, and are the eigenvalues of B^, a, gj > gj, and Kcep„Bj is a function 
of (1-Pmd) probability and [o^^sj/^i^Bjj^^ in the CEP look-up table. 

The total estimated horizontal position error for the jth Sub-Filter is computed as 

HPEj = HPE_Bj + HPE^NPj - HPE_NBj (44) 
Since there are multiple Sub-Filters, HPL^i is determined as follows: 

HPLHi=max HPEi (45) 

HPL Based on the Rare Normal Fault-Free Hypothesis 

If there is no fault in the navigation solution, the Kalman error covariance from the 
Full-Filter can be used to determine HPL^o by following steps (1) and (2). 

HPLho=K^.V^ (46) 

where a, ^ and 03 0 are the eigenvalues of , a^ Q> o^ ^ , Kffd is a function of (1- 

Pho) probability and (cr, 0/^2,0/'^ *he CEP look-up table, and Pho is the fault-free integrity 
probability. 
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HORIZONTAL UNCERTAINTY LEVEL 

The algorithm for the Horizontal Uncertainty Level is adapted from the algorithm 
developed for the HPL. Therefore, like the HPL, HUL consists of two components— HUL^o 

and HULhi . By proceeding with the same argument in developing Eq. (46), it can be 

shown HULho = HPLjjo . 

Horizontal Uncertainty Error Due to Bias and Uncorrelated Noise 

Since X-(k) contains the information about the 6iaso(k) vector, it can be used to 

estimate the portion of the true horizontal position error due to the bias and uncorrelated 
noise at time k. By proceeding with arguments similar to those used for Eq. (34), the 
horizontal uncertainty error due to the bias^ik) vector can be computed for a specified 
probability as follows: 

HUE_B^(k) = ||AiaSo(k)|| 

(47) 

where Apj is intended to compensate for the reduction of y^(k) by noise in the y 
domain. 

Now. the question is how to determine Apj such that Eq. (47) will bound the current 

horizontal position error induced by the bias and uncorrelated noise with the (1-Pmd) 
confidence. As derived in Eq. (28) of Ref. [5], Apj can be approximated as follows: 

Apj « Pbias-Td (48) 
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where Pbias is given in Eq. (32) and Td is the detection threshold used in the fault 
detection. 

So, Eq. (47) can be rearranged as follows: 

HUE^B.(k) < V<^,Bj(k)(VMk) + Pbias-Td) (49) 

Horizontal Uncertainty Error Due to Noise Only 

The horizontal position errors due to noise derived in Eqs. (42) and (43) can be used 
to compensate HUE_Bj(k). Therefore, the total estimated horizontal uncertainty error for 

the jth Sub-Filter can be computed as follows: 

HUEj(k) = HUE_B-(k) + HPE_NPj(k)-HPE__NBj(k) (50) 

Since there are multiple Sub-Filters, the HULHi(k) is determined as follows: 



HULhi (k) = max HUE^ (k) (5 1 ) 

i=l 



Therefore, HUL at time step k is determined as follows: 

HUL(k) = max{HULjjo(k), HULHi(k)} (52) 

FAULT EXCLUSION PROCESSING 

In this embodiment of the present invention, measurement residual monitoring is 
proposed for the fault exclusion function. There are two requirements for fault exclusion. 
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(I) Fault identification: Identify the failed measurement with the probability of (1- 
Pmi), where Pmi is the specified probability of misidentiflcation. 

(II) Post- 
exclusion fault detection: If the failed measurement is excluded, the fault detection 
function based on remaining measurements should be available for the given phase of 
flight. 

Now referring to Figure 5, the basic processing flow for the proposed fault exclusion 
function is shown. 

The blocks with bold solid lines are designed to meet the first requirement for fault 
exclusion, and the block with bold dashed lines is designed to satisfy the second 
requirement. The function of each of these blocks in Figure 5 will be outlined below. 

Fault Identification 

As shown in Figure 5, there are two independent fault identification processes — 
post-update residual monitoring and least-squares fault detection and identification. The 
former is primarily designed to isolate medium-to-fast failure rates, and the latter is used to 
isolate slow failure rates when there are six or more measurements available. Since these 
two processes are complementary to each other, they are combined here to provide better 
availability for the fault exclusion function. Either of these two processes can trigger the 
correct fault isolation and move forward to the next step. 



-28- 
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Post-Update Residual Monitoring 

The fault isolation algorithm presented in Parkinson, B. W., Axelrad, P., 
"Autonomous GPS Integrity Monitoring Using the Pseudorange Residual", NAVIGATION, 
Journal of the Institute of Navigation, Vol. 35, No. 2, Summer 1988, pp. 255-274 hereafter 
referred to as Ref. [10] is adapted below for our applications. The post-update residual 
vector at time k for each Sub-Filter is computed as follows: 

^(k) = z,ik)-H,(k)x, (k/k), i=l,...,N (53) 

where is the post-update residual vector for the ith Sub-Filter and is the 
linearized measurement vector for the ith Sub-Filter. 

For applications related to Kalman filters, the innovation, z^{k), is generally used to 
exclude abnormal measurements Ref. [3], where 

zr(k) = z,(k)-^,(k)x,(k/k-l) (54) 

However, for the following reasons, the post-update residual vectors are used 
instead. 

• The post-update residual is a white sequence and its covariance is smaller 
than the covariance of the innovation. 

• Due to the growth of the clock bias uncertainty between Kalman filter 
updates, the innovation variance can be relatively large. Therefore, it is more difficult to 
detect the presence of failures with the innovation than with the post-update residual. 
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• The primary advantage of the innovation over the post-update residual is that, 
for pseudorange step faults with magnitude more than several hundred meters, the fault 
can be detected and excluded before the failed measurement is incorporated into 
solutions. For this reason, the step detector function found in most aviation GPS 
receivers uses innovation to exclude large step errors. However, for ramp failures and 
smaller step errors, the step detector function will not respond and the FDE function 
must be relied upon. For these kinds of failures, the FDE function would generally 
detect the presence of abnonnal solutions after failed measurements have been 
incorporated into solutions. Therefore, there is no real advantage in using innovation 
over the post-update residual, especially for excluding slow ramp failures. 

The post-update residual vector for the jth Sub-Filter not using the failed 
measurement can be shown to be zero mean and have a covariance, (k) , as given by 



4(k) = if.(k)-^.(k)P.(k/k)fl^/(k) 



(55) 



Now, we define the residual test statistic, 7j (k), for the jth Sub-Filter as 



y,(k) = z7(k)4-(k)^.(k) 



(56) 



The direct computation of A^^ (k) in Eq. (56) can be avoided by computing the 
Cholesky decomposition, 

4(k) = Z..(k)£j(k) (57) 
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where Xj (k) is a lower-triangular matrix. The residual test statistic is then computed 

as 

C,(k) = L:>(k)z.(k) (58) 

Yi(k) = C/(k)C3(k) (59) 
where Eq. (58) is computed via forward substitution. 

It can be shown that Yj(k) is central chi-square distributed with (N-1) degrees of 
freedom since the jth Sub-Filter does not use the failed measurement. It can be further 
shown that (k) for i ^ j would be noncentral chi-square distributed with (N-1) degrees of 
freedom because all other Sub-Filters use the failed measurement. This difference in the 
statistical distribution between yj (k) and yi(k) is used as the basis for the fault 

identification. The required correct fault identification with the probability of (I-Pmi) for a 
particular j is 

(y.y''<T,and(y,fST„ j, i=l,...,N (60) 

where \ is the isolation threshold and the determination of (Tj)^ is based on the 
central chi-square statistics with (N-1) degrees of freedom and probability of (1-Pr^4i), is 
the exclusion threshold whose determination is discussed later. 

Since the post-update residual vector can be shown to be uncorrelated with 
previous post-update residual vectors, the performance of the fault identification process 

31 - 
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can be further improved by averaging the post-update residual vector over several 
measurement update cycles to reduce noise. However, as mentioned in Ref. [8], the 
averaging technique would result in the reduction of the detection performance for failures 
with dynamics that are fast relative to the averaging period. Therefore, to isolate failures 
due to different failure rates, the post-update residuals are averaged over three different 
intervals (in terms of update cycles) — 1 , 10, and 30. So, there are three averaged residual 
parameters in all for each Sub-Filter. 

As shown in Figure 5, the fault detection function and the fault exclusion function are 
two independent processes that run in parallel. For fast failure rates, this approach would 
enable the post-update residual monitoring method to exclude the failed measurement 
even before a fault is detected in the position domain because the large transient errors 
induced by fast failure rates would show up earlier in the range domain. However, with the 
fault detection and the fault exclusion as two independent processes, we need to ensure 
the fault exclusion processing can meet the same false alarm requirement as the fault 
detection processing does. Because there are three averaged residual parameters for 
each Sub-Filter and any of them can cause a false alert, needs to be determined with 

the false alert probability of Pfa/3. Therefore, the determination of (T^f is based on the 

central chi-square statistics with (N-1 ) degrees of freedom and probability of (1 -Pfa/3). 
With Tj and determined in this manner, Eq. (60) can ensure the jth Sub-Filter is fault- 
free with (1-Pmi) probability, all other Sub-Filters are non-central chi-square distributed, and 
the false alert requirement is met. 
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Least-Squares Fault Detection and Identification (LS FDE) 

It has been found that there is one drawback in using the post-update residual 
monitoring to perform fault identification. Even with the averaging technique, this form of 
monitoring may not be able to isolate failures induced by slow pseudorange failure rates 
(e.g., 0.02 m/s or lower). However, when there are six or more measurements available, 
the following LS FDE method can be used to complement the post-update residual 
monitoring method. 

As shown in Young, R. S. Y., "Oceanic/Remote FDE Algorithm Design", internal 
memo, Rockwell Collins Inc., August 16, 1996 and Lee, Y., "New Techniques Relating 
Fault Detection and Exclusion Performance to GPS Primary Means Integrity 
Requirements", ION GPS-1995, Palm Springs, CA, Sept. 12-15, 1995, pp. 1929-1939 
hereafter respectively referred to as Refs. [1 1] and [12], when there are at least six 
measurements available, the snapshot algorithms can be used to perform fault isolation. 
Since the performance of the snapshot algorithms is independent of failure rates, they 
effectively complement the post-update residual monitoring method. The fault isolation 
algorithm in Ref. [1 1] requires less computational time since they do not require multiple 
navigation solutions. 

The drawback of using snapshot methods for isolating slow ramp failures is that it 
requires at least six measurements. Therefore, when there are less than six 
measurements available and the pseudorange failure rate is slow, the fault identification 
process may not be able to isolate the failed measurement. In this situation, the HUL 
should be used to monitor the horizontal position error. If HUL is greater than the HAL for 
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the given phase of flight, the navigation solution from the integrated GPS/inertia! sensor 
should not be used. OthenA/lse, the navigation solution is tolerable and can still be used. 

Post-Exclusion Fault Detection 

The second requirement for the fault exclusion indicates that, after the failed 
measurement is excluded, the HPL based on the remaining measurements (post-HPL) 
should be less than the HAL for the given phase of flight. The primary purpose of this 
requirement is to ensure the integrity function can detect another measurement failure that 
occurs after the first failure is isolated. 

If the failed measurement is correctly identified and excluded, the jth Sub-Filter 
omitting the failed measurement should be fault-free. Therefore, the solution from the jth 
Sub-Filter can be used to re-initialize the Full-Filter and other Sub-Filters. However, before 
the failed measurement can be excluded, the second requirement for fault exclusion needs 
to be met. At the time of fault isolation, the post-exclusion HPL is the HPLho from the jth 
Sub-Filter since all filters will be initialized with the same solution. So, if the HPLho from the 
jth Sub-Filter is less than the HAL for the given phase of flight, the failed measurement can 
be excluded. 

After the re-initialization, the Full-Filter would only use (N-1) measurements and 
each Sub-Filter would only use (N-2) measurements. The post-HPL is computed from the 
HPL algorithm proposed before. 



PATENT 

Attorney Docket No. 02CR109/KE 



A more thorough understanding of the present invention may be achieved by 
reviewing the more detailed description below regarding the covariance of the solution 
separation vector, (k) . 



The computation for the covariance of the solution separation vector, *j (k) , defined 

in Eq. (8) will be derived here. Since the jth Sub-Filter uses one less measurement than 
the Full-Filter, there is one pseudorange bias error state used by the Full-Filter, but not by 
the jth Sub-Filter. However, to simplify the derivation, it is assumed the jth Sub-Filter still 
carries this unused pseudorange bias error state. It can be shown later that carrying an 
unused pseudorange bias error state, which is modeled as a first-order Markov process, 
does not affect the computation of (k) and the performance of the jth Sub-Filter. 

The time evolution of estimation enrors for the Full-Filter and jth Sub-Filter is shown 

as follows: 



When the filter processes measurements, these two errors are changed according 

to 



COVARIANCE OF THE SOLUTION SEPARATION VECTOR 



Jco(k/k-l) = 0(k.k-l)xo(k-l/k-l)-(D(k) 



(A1) 



X. (k/k-1) = C>(k,k-l)Xj (k-l/k-1) - <o(k) 



(A2) 



io(k/k)=(/-j!ro(k)^o(k))xo(k/k-i)+jj:o(k)uo(k) 



(A3) 
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Xj(k/k) = (/-jr.(k)^.(k))jc.(k/k-l) + J5:j(k)u.(k) (A4) 

where JSToCk) and K^Qa) are the Kalman gains for the Full-Filter and the jth Sub- 
Filter, respectively. 

By combining Eqs. (A1)-(A4), the evolution of the cross-covariance Poj(k/k) can be 
computed as follows: 



Po,.(k/k) = To(k,k-l)P,^.(k-l/k-l)^^(k,k-l) 



(A5) 



where 

%ik,k-l) = [I- JS:o(k)^o(k)]<I>(k,k-l) (A6) 

^"(k,k-l) = 0^(k,k-l)[/-ii:j(k)^j(k)]' (A7) 

r,ik) = [l-K,(k)H,ik)] (A8) 

r](k) = [l-K.(k)H.ik)y (A9) 

i?^.(k) = E{uo(k)oJ(k)} (A10) 

The observation matrix of the jth Sub-Filter, (k) , is related to the observation 
matrix of the Full-Filter, Hf,{k), as follows: 

H,{k) = E,H,{k) (A11) 

-36- 
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The measurement covariance of the jth Sub-Filter, (k) , is related to the 
measurement covariance of the Full-Filter, (k) , as follows: 



With Eq. (A12), we can see that the jth column and jth row of Rj{k) are all zeros. 



where if/(k) denotes the Moore-Penrose generalized inverse of the ifj(k) matrix. 

With Eq. (A13), we can see that the jth measurement is not used by the jth Sub-Filter. 
Therefore, carrying the pseudorange bias error state for the jth measurement has no effect 
on the performance of the jth Sub-Filter. 

The initial value for can be determined as 



*i(k) = ^A(k)^j 



(A12) 



K.(k) is computed as follows: 



jr,.(k) = i>,(k/k)^/(k)j?;(k) 



(A13) 



(0/0) = E {(xo (0/0) -x{0)){x. (0/0) - X {0)f } 
= E{x, (0/0)xJ (0/0)) -E{:c(0):c^ (O)} 



(A14) 



With the Initialization specified by Eq. (3), Eq. (A14) can be re-arranged as 



Poj (0/0) = E {io (0/0) xl (0/0)} -E{x{0)x' (0)} = P, (O/O) 



(A15) 



With Eqs. (4) and (A15), we can show 
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Poj(0/0) = i> (0/0) = Pj(0/0) = Po (A16) 

From Kalman filter theory, we know 



i»(k/k)=[/-jro(k)^o(k)]* 

[o(k,k-l)Po(k-l/k-l)<I>^ (k,k-l) + e(k)] 



(A17) 



Define 

A,{k) = K,{k)H,{k) (A18) 
Witli Eqs. (A16)-(A18), Eq. (A5) can be re-arranged as follows for k = 1. 
Pp, (1/1) = Po(l/l) - Po(l/lK (1) 

JSTq (l) and (l) in Eq. (A1 9) can be shown as: 

K,{l) = P,{in)Hl{l)R;\l) (A20) 
S,,{l) = S,{l)E, (A21) 

Now, let us examine the last two terms at the right hand side of Eq. (A19). 

= [-Po(l/l) + Po(l/l)]4(l) (A22) 
= 0 

With Eq. (A22), Eq. (A19) can be simplified as: 

-38- 
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P^(l/l) = Po(l/l) (A23) 

Now, the process between Eq. (A19) and Eq. (A23) can be repeated for k = 2,3... 
Therefore, by induction, Eq. (A23) is valid for all k. 

With Ppj(k/k) = Po(k/k) established for all k, ^^(k) can be computed as 

J»- (k) = (k/k)- (k/k), k = 0,U2„.. (A24) 

where Pj^ (k/k) is the matrix containing the first 2x2 elements of Pj (k/k) and 
Pq (k/k) is the matrix containing the first 2x2 elements of P^ (k/k) . 

In Ref. [4], a similar result as Eq, (A24) was derived for the more restricted case — a 
Full-Filter using all measurements and a Sub-Filter using no measurement. Since /^-^ (k/k) 

and Pq (k/k) are already available from the normal Kalman filter processing, (k) can be 

easily obtained as the difference of these two covariance matrices. Therefore, the sf)ecial 
dual covariance propagator used in Refs. [1] and [2] to compute is not required. 

Due to Eqs. (A16) and (A24), there is a singular point (6(0) = 0) for the fault 
detection function at initialization. This characteristic is expected since there is no 
difference between Xq(0/0) and jCj(0/0). 

A more thorough understanding of the present invention can be achieved by 
reviewing the following conclusory statements: A real-time fault detection and exclusion 
algorithm is presented for tightly integrated GPS/inertial sensors. The fault detection 
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algorithm is based on normalized solution separation. With normalization, the proper 
probability levels for the detection threshold and HPL can be obtained. In addition, HPL is 
predictable without real GPS measurements, which is an important operational 
requirement. 

A hybrid algorithm that combines two independent methods that complement each 
other is proposed for the fault exclusion function. For medium-to-fast failure rates, the 
post-update residual monitoring scheme performs better and does not require a minimum 
of six measurements. For slow failure rates, the least-squares fault identification scheme 
performs better when there are at least six measurements available. In addition, since it 
takes an extremely long failure exposure time for slow ramp failures to corrupt the 
navigation solution to be above HAL, the continuous constellation change due to satellite 
motion may provide enough observability during some periods of the long failure exposure 
time to permit the correct fault isolation. During the periods when the correct fault isolation 
cannot be performed, HUL can be used to monitor the horizontal position error. 

It is thought that the system, apparatus and method of the present invention can be 
understood from the above description. It should be understood that the above description 
is not exhaustive and is an example of one of many variations of the present invention 
which could be performed based upon the spirit and scope of the invention. It is the 
intention of the inventors to claim all such variations within the appended claims. 



